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We present an algorithm for the detection of periodic sources of gravitational waves with inter- 
ferometric detectors that is based on a special symmetry of the problem; the contributions to the 
phase modulation of the signal from the earth rotation are exactly equal and opposite at any two 
instants of time separated by half a sidereal day; the corresponding is true for the contributions from 
the earth orbital motion for half a sidereal year, assuming a circular orbit. The addition of phases 
through multiplications of the shifted time series gives a demodulated signal; specific attention is 
given to the reduction of noise mixing resulting from these multiplications. We discuss the statistics 
' of this algorithm for all-sky searches (which include a parameterization of the source spin-down), 

I in particular its optimal sensitivity as a function of required computational power. Two specific 

■ examples of all-sky searches (broad-band and narrow-band) are explored numerically, and their per- 
04 ' formances are compared with the stack-slide technique (P. R. Brady, T. Creighton, Phys. Rev. D 

\ 61, 082001). 
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^ ; INTRODUCTION 

a^: 

^-H , Kilometer-scale gravitational wave interferometric detectors will become operational by the end of 2001, allowing 
' observations with unprecedented sensitivity at frequencies ranging from a few tens of hertz to approximately one 
, kilohertz ||l|. Periodic sources form an interesting subclass of potentially detectable astrophysical objects, due to the 
possibility to improve the "visibility" of a weak periodic signal in noisy data by increasing the observation time. Good 
candidates for such sources are spinning galactic neutron stars with a deformation misaligned with their rotation axis; 
optimistic estimates suggest that LIGO I could detect such objects with an observation time of the order of one year 



Frequency modulations due to the Doppler shift induced by the detector motions, and intrinsic frequency variations 
due to the loss of angular momentum to gravitational waves, are two easily parameterizable examples of properties 
of a realistic signal that will dramatically increase the difficulty of the data analysis. In fact, for all but the simplest 
^ ■ types of searches, the data analysis will be sub-optimal because of computational limitations. More complications 
k>( \ are likely to arise if the source is in a binary, due to the additional frequency modulation and possibly larger or 
j_j i random intrinsic frequency variations expected if the source is accreting. 

' These considerations have driven serious efforts to define data analysis schemes that would maximize the sensitivity 



of a search (to be defined in section IB), given a maximal available computational power. Matched filtering analysis, 
together with an in-depth discussion of signal properties, is presented in a series of papers [Q; for all-sky searches, this 
technique would be computationally prohibitive. The so-called coherent search, discussed in uses a resampling 
technique to demodulate the signal and to correct for intrinsic frequency variations: the detection is performed by 
taking the Fourier transform with respect to the resampling time, which amounts to a change of frame of reference 
from the detector to the solar system barycenter and to a "stretching" of time for frequency variations. This technique 
is generally not practical for all-sky searches over long time periods (longer than a few days), but a refinement of 
it called the stack-slide technique [g| is considered as one of the best of all known techniques: the initial data are 
segmented into shorter intervals, and their power spectra are built using the coherent method. For these shorter 
intervals, a smaller number of points in parameter space have to be used. The spectra are added (incoherently) after 
having been shifted according to a finer gridding of parameter space. A related technique uses Hough transforms 
to track the frequency evolution of peaks in the spectra from the shorter intervals, which is known in advance for some 
choice of the source parameters. In general, it is possible to improve the sensitivity of these techniques by integrating 
them into more general hierarchical searches : a first search in parameter space is done at reduced sensitivity but 
low confidence level (so that spurious events are hkely); the vicinity of candidate events is then scrutinized at full 
sensitivity to get higher confidence detections. 

We have explored yet another approach that exploits a particular symmetry of the problem in an attempt to reduce 
its complexity. Firstly, we note that the Doppler shift in the frequency of a given source due to the rotation of the 



1 



earth, independently of the source position, is exactly equal and opposite at any two instants of time separated by 
half a sidereal day (we will only consider linear terms in the Doppler shift, an excellent approximation). Secondly, 
in the approximation that the earth orbital motion around the sun is circular, the Doppler shift in frequency has 
the same property as above, but for a time separation of half a sidereal year. We shall use this approximation for 
illustrative purposes; corrections for the small eccentricity of the earth orbit and for the influence of the moon and 
of other planets would have to be included in practice. Signal phases at different times can be added by properly 
multiplying the (analytic extensions of the) signal time series, and therefore a signal free of frequency modulations 
may in principle be built. We will give a more precise formulation of this idea in section |. 

The multiplications of different stretches of the data will strongly affect the noise characteristics, principally by 
mixing noise at different frequencies; a useful implementation will have to bandpass filter the data to minimize 
this effect, and therefore a search over source parameters (frequency, frequenc y de rivatives, source position) will be 
required. An example of a possible implementation will be presented in section 1 A , The achievable sensitivity (to be 
defined in section [ B ) will be studied numerically as a function of computational power requirements in section ||, 
and comparisons with the stack-slide method mentioned above will be presented. 



I. ALGORITHM 



We will use the following model for the signal s{t) at the detector: 



s{t) — hcos[(j){t)], 
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where e is the earth oblicity, A is the detector latitude, c is the speed of light, i?® is the earth radius, A.U. is an 
astronomical unit, d is a sidereal day length, yr is a sidereal year length, (a, S) is the source angular position, /o is the 
signal frequency in the source rest frame (it may be a function of time), and is a fixed number that defines the origin 
of time. For simplicity, since our focus is on the phase modulation of the signal, the explicit time dependence of h on 
the detector response pattern will not be considered, although its averaged effect will be included in our definition of 
sensitivity (section ^^). If /o is independent of time, a good approximation to the signal half-bandwidth A/ over a 
time interval of length t (i.e., half the intrinsinc frequency change over t), is given by: 

A/(0 =/o[Aimin(l,</d)-}-A2min(l,t/yr)]. (7) 

The numerical values of the terms multiplying the cosine terms in eq. ^ and eq. ^ are approximately 1.55 • 10"^ and 
9.94 • 10~^, respectively. Hence, the half-bandwidth of a 1 kHz source would be roughly 2 mHz over one day, and 0.1 
Hz over a year. 

The phase function (p{t) is sufficiently well-behaved that we can use the Hilbert transform of s{t) to build its 
quadrature s{t) — /i sin[(/)(t)], and therefore construct the analytic signal S{t): 

S{t) = he"^^'\ (8) 
Assuming that /o is independent of time, we have that 

Si{t) = S{t)S{t + d/2)5(t + yr/2)S{t + d/2 + yr/2) 

= /j4g87ri/o(t+d/4+yr/4)^ (9) 

that is, ^4 is monochromatic with frequency 4/o. If /o is a function of time, however, 5*4 as defined above will not be 
demodulated because of the coupling of /o with the modulations; its phase will be 
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^ - fo(t)t + + k{t2)t2 + k{h)h + 

ZTT 

Ml{t)[fo{t) - /o(tl)] + Afi(i2)[/o(t2) - h{h)] + 

M2{t)[h{t) - /o(i2)] + A/2(ti)[/o(ti) - /o(i3)], (10) 

where = t + d/2, ^2 — ^ + yi'/2 and = t + d/2 + yr/2. The signal so obtained stiU presents some residual phase 
modulation from the detector motions, but it results in a much smaller contribution A/mod to the signal bandwidth; 
the terms multiplying the modulation functions Mi and M2 are now of order of the intrinsic frequency change over 
half a day and half a year, respectively, rather than the source frequency itself. We will show below that our algorithm 
best performances are obtained by dividing the total data from an observation of length iobs into smaller segments 
of length T shorter than one day. Consequently, we impose the condition that the signal bandwidth produced by the 
residual modulations in eq. |l^ be smaller than the frequency resolution for a segment of length T, A/mod < 
assuming a linear frequency model: /o(i) = /o • [1 — ^/t], where r is the source spin-down characteristic time. We find 
that this condition is: 

T> kT\Ai+A2), (11) 



^> (tI^V^ ) -SSyr. (12) 
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This estimate shows clearly that for most realistic sources, the demodulation by multiplication is performed correctly, 
so that we can assume that spin down enters the phase function only through the combination fo{t)t + fo{ti)ti + 

k{t2)t2 + k{h)h. 

Noise will be seriously affected by the multiplications of data segments at different times, through nonlinear mixing 
of signal with noise and of noise with noise at different frequencies. The latter of these effects is so severe that we will 
have to minimize it as much as possible by bandpass filtering individually the data segments around the signal, or the 
technique will be of no utility. An efficient filtering requires a knowledge of the frequency content of the signal as a 
function of time, and will therefore involve a search over source frequencies, sky positions and spin-down parameters; 
this will greatly increase the computational costs. This search will be implemented by meshing this parameter space 
with a certain resolution, i.e. with a certain number of points Np. Obviously, if a given signal has parameters that 
are not exactly coincident with any of these points, the sensitivity of our search to this signal will be reduced. We will 
explore below how the density of the mesh couples to the achievable mean sensitivity, and therefore give an estimate 
of the usefulness of our technique. 



A. Implementation 

From the general guidelines established in the previous section, we are able to discuss an implementation of our 
algorithm. This implementation has parameters that allow an optimization of the method, subject to computational 
power constraints. 

The first step will be to prepare for the analysis data obtained from an observation of length ^obs- The dataset 
is compressed to contain only information for frequencies of interest, i.e., frequencies between /min and /max- The 
analytic signal is then constructed from this reduced dataset, using Hilbert transforms to build the quadrature, as 
described in section |. The resulting data is divided into equal segments of length T. The computational cost of these 
operations is negligible, especially because they are only performed once at the beginning of the search. 

Next, a bank of bandpass filtered segments is built. We choose a filter half-bandwidth B > A/, where A/ is the 
maximum of the signal half-bandwidth over time T, and for every segment, we construct Nf smaller filtered segments 
by using a heterodyne technique pO[ . These filtered segments are arranged to cover completely the frequency interval 
from /min to /max- If wc imposc the requirement that no signal power should be lost due to the filtering, the number 
Nf is a function of -B — A/: the whole frequency interval is divided into subintervals of length 2B, displaced from 
each other by 2{B ~ A/), so that any signal with half-bandwidth A/ is guaranteed to be fully contained within one 
subinterval, and Nf = (/max — /min)/2(i? — A/). In the heterodyne technique, the signal is multiphed by a complex 
exponential exp 27ri(_B — /o)t, low-pass filtered (for reasons of antialiasing), and resampled at frequency 8B. The 
resulting segments have 8BT (complex) data points, and their Nyquist frequency is 4B. The cost of the filtering, if 
implemented by an infinite impulse response digital filter, is of order 2(2n — 1) fioating point operation per input point, 
where n is the number of poles in the filter, which is a function of the filter performance; we shall denote this cost by 
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Cf. Resampling will have negligible computational cost compared to the multiplication by the complex exponential; 
assuming that the latter is not computed at run time, the cost for the construction of the filtered segments bank is 

I /~i 1 c\-i (/max /mill) n o^ 

(Cf + 6)tobs 2{B-Af) 

This cost can be adjusted by varying B and T, which enters the formula through the dependence of A/ on it. 

Once the filtered segments bank is constructed, for each of the Np points in parameter space, and for every segment 
at time t < tobs — T — d/2 — yr/2, we select the appropriate filtered segments from the bank, i.e. filtered segments 
where the signal is expected for times t, ti = t -\- d/2, t2 = t + yr/2 and = t + d/2 + yr/2, and we multiply them 
together. At this point, the modulations induced from the detector motion have been removed; what remain to be 
applied are corrections for the spin-down of the source. A refined gridding of the spin-down portion of parameter 
space, for the already quite narrow region used for the choice of filtered segments from the bank, is used to resample 
the product of the four segments according to the new time 

/o(0) ' ^'^^ 

where the origin of t was arbitrarily set to zero. The resampling will only affect a few points in the time series, so its 
cost is negligible in itself. The power spectrum is constructed by taking the norm of the Fourier transform. Finally, 
all the spectra corresponding to a given point in parameter space (from the partitioning of time in segments of length 
T) are added together to form the final spectrum, which is searched for significant peaks up to the Nyquist frequency 
AB. All these operations are done on small segments of length 8BT, and the total cost is the sum of the cost for 
building the products for all Np points and of the cost for the spin-down corrections, power spectra construction and 
sum for the Np points corresponding to the finer grid: 

Btob^Np[72 + 2QN'p{l + Xog^ 8BT)] flop. (15) 

Again, this can be adjusted by varying B and T, since Np and Np depend on both B and T. The total cost of the 
search is the sum of equations ^ and 

The number Np will be the number of non-intersecting (hyper-) volume elements in parameter space required to 
cover completely the region of interest, such that all points within a given volume element correspond to the same 
choice of filtered segments from the bank, at each of the four times of interest in this problem. Manifestly, Np will 
depend on B through the way it partitions the frequency interval [/min, /max], and on T, which influences the size of 
a volume element over which the source frequency varies by some flxed amount. We solve the problem of computing 
Np numerically, by building a mean volume element: we select at random a point in parameter space, and compute 
its associated volume element. Repeating this procedure, we construct the average of a volume element, and divide 
the parameter space volume by it; the result is an estimate of Np. 

The refined gridding should only be necessary for searches with large spin-down rates (see below), and its effect 
will be to force to repeat Np times the operations following and including the resampling above, where Np is the 
number of points on the refined grid. The computation of A^^ is not of the same nature as the one for Np-. the 
problem now only involves the spin-down parameters, and the goal is to minimize the losses in signal power from 
residual frequency drifts, or misalignments of the spectra that are summed together. In that sense, it is similar to 
the problem of computing the number of points in parameter space for the stack-slide technique and should be 
solved using the same geometrical approach. However, we shall rather take Np = 1, and justify this approximation 
by the relatively large spin-down time (r = 5000yr) to be used in the numerical example of section It should be 
kept in mind that this may not be sufficient for faster spin-down, in the sense that under this approximation Np may 
be under-estimated. With the ability we now have to compute the required computational cost as a function of B 
and T, the only thing that remains to be described is how we characterize the sensitivity of the search, as a function 
of the same choice of parameters. 



B. Sensitivity 

From a statistical point of view, the results in |^,^,^ are presented in a somewhat non-standard way, so to ease 
comparison of our results with theirs, we first briefiy reformulate their definition of sensitivity. By setting the average 
noise power to be equal to one, we can express the signal amplitude h in units of yo^XfV^jbs, where S'„(/) is the 
noise spectral density at frequency /. In order to be consistent with |§,||J^, we use the definition that the averaged 
signal power (also the square of the signal to noise ratio in amplitude) is Ps = {F^)h'^, where (F^) = 1/5 is the average 
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of the detector response over all angles and polarizations This definition is useful to approximately account for 
variations in h over time that are produced by the motion of the source in the detector response pattern. For a true 
signal power Pg, computing the power spectrum of a time series as a sum of n of its subseries spectra (this describes 
the statistics of both the coherent and stack-slide techniques) will give a distribution pn{P\Ps) of the power P 
Given a false alarm probability 1 — a, one defines a threshold power Pa such that the integral of»„(P|0) from to Pa 
is a. Observations with power exceeding Pa will be called "detections with confidence level a" Q. One further defines 
a false dismissal probability /3 as a function of the signal power, as the integral of pn{P\Ps) with respect to P, from 
to Pa ■ It has been customary in the field [§|j5|J^ to use the relation (P) — Pg +n for the average observed power (P) 
to interpret Pa as a lower bound on the detectable signal amplitude, by setting Pa — (P). For instance, performances 
of sub-optimal techniques are frequently compared to the optimal case defined by n = 1; for a = 0.99, the minimal 
detectable amplitude as defined above (the optimal amplitude) is h — 4.2, and the corresponding false dismissal rate 
is /? ~ 0.57. In fact, it is easily checked that as n gets larger, this definition gives a value of /3 that approaches 1/2 
from above. Therefore, for the ease of comparison with other works, we will define our minimal detectable amplitude 
as the smallest value of h that corresponds to a false alarm probability 1 — a and a false dismissal rate /? — 1/2. The 
sensitivity will be defined as 8 = 4.2/ft,, and this definition will be approximately compatible with those of 

We use numerical analysis to determine ft, as a function of B and T. We first set the signal amplitude to zero, 
and generate a large number of filtered segments of pure white noise, that are multiplied together and then added by 
groups of (tobs — T — d/2 — yr/2)/T (i.e., the number of segments in the portion of the observation period that can 
be used for multiplications). By constructing the cumulative probability function for the resulting segments, we then 
evaluate the threshold Pa- We repeat the procedure for non-zero signal amplitudes: picking at random in a uniform 
distribution the parameters of the source, we construct the final spectrum, and evaluate the power in the expected 
frequency channel. For each choice of signal amplitude, we repeat this a large number of times for different source 
parameters, and build a distribution of the power in the expected channels; the values of /9 as a function of signal 
power are then deduced from these distributions and the previously calculated value of Pa- Therefore, the sensitivity 
O is effectively averaged over all possible source parameters. To reduce computational burden, spin-down effects are 
not directly included, in the sense that we only consider them in determining the bandwidth of the signal for a given 
T, and do not include them in the initial simulated segments that are multiplied together. By doing so, we might 
reduce the non-linear mixing of signal with noise, and therefore deduce a sensitivity that is too high; however, the 
magnitude of this mixing is such that the error so induced is expected to be small. 

II. NUMERICAL SIMULATIONS 

We opted for a relatively "easy" broad-ha,nd all-sky search for our explorations. We choose /min = 40Hz, /max = 
IkHz, T — 5000yr, iobs = lyr, and confidence level a = 99%. Note that iobs = lyr is the minimal possible observation 
time, and that the data at times t > d/2 -|-yr/2 can not be demodulated. This obviously limits the sensitivity of the 
search, although it does so in a diminushing proportion as iobs grows above one year. A different narrow-band all-sky 
search was also considered, with parameters as for the broad-band search, except for = 450IIz and /max = 500IIz, 
perhaps corresponding to the case of an interferometer made narrow-band using signal recycling. For definitiveness, 
we choose the latitude of the detector to correspond to the LIGO Hanford detector. As described above, we computed 
numerically the value of Np for different values of T and B (from now on, B will be expressed in units of A/(T) as 
defined in eq. 0, with /o = /max, and Ai, averaged over all angles 5); see figure |l|. Using the method described in 



section [B, we computed the dependence of the relative sensitivity, O, on the same variables; the results are shown 
in figure]^ 

To get a meaningful notion of the computational power P from the cost of the search (the sum of equations |l^ and 
p^ ), we imposed the requirement that the analysis should be done in a time iobs- Our map of TVp vs. B and T was then 
reexpressed as a map of P versus B and T. We observed that for plausible filter implementations, the contribution to 
P from eq. |l^ was at least two orders of magnitude larger than the contribution from eq. |l^ in the region of interest of 
the r, i3-plane; our results are highly insensitive to the real cost of constructing the filtered segments bank. Finally, 
we optimized our algorithm for a given computational power by choosing the values of B and T that maximize O 
along a line of constant P. The results are shown in figure ^, together with the corresponding sensitivity for the 
stack-slide technique. 



^The authors of j^Q define the false alarm probability as 1 — ct/Np, 
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FIG. 2. A contour plot of S as a function of B (in units of eq. M) and T. Labels on the contour lines are values of Q. 
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FIG. 3. Top: the relative sensitivity O for the optimized algorithm, as a function of computational power P (solid line), and 
the corresponding for the stack-slide technique (dash-dotted line), both for the broad- band search. The dotted line corresponds 
to the different narrow-band search. Middle: the optimal segment length T in the broad-band search. Bottom: the optimal 
filter bandwidth B (in units of eq. also in the broad-band search. 
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A. Discussion 



Fig. H shows that the sensitivity of the algorithm described in this paper is roughly 30% of the sensitivity of the 
stack-slide technique in a broad-band search, for computational power ranging from 0.4 Tflops to 8 Tflops. For the 
same range of computational power, the minimum amplitude detectable with 50% efficiency ranges approximately 
from 25 to 60 times the optimal amplitude. 

Also shown in fig. ^ is the sensitivity of an easier narrow-band all-sky search. When the computational power 
reaches 4 Tflops, the sensitivity eventually gets better than that of the stack-slide technique, as computed for the 
harder, broad-band search. However, it is not clear how the performances of the stack-slide technique scale with the 
bandwidth of the search, although they do scale strongly with the upper frequency. It is not possible, for instance, to 
apply the stack-slide technique on a heterodyned dataset without having to search over the additional source frequency 
parameter, because of the coupling of the source frequency to both the detector rest frame time and to the phases 
introduced by the detector motion (Mi and M2 in equation ^). The comparison of our algorithm to the stack-slide 
technique is therefore not trivial, and in particular there might be some narrow-band searches where the stack-slide 
performances are approached or outperformed. 

Our results suggest that the distortions in the noise produced by the nonlinearities inherent to our algorithm can 
not be alleviated by narrow-banding the signal such that the algorithm performs in a manner that is satisfactory for 
a realistic broad-band all-sky search. The need for narrow-banding adds one dimension (the source frequency) to the 
parameter space, and this greatly increases the computational burden. Moreover, the dependence of the sensitivity 
on the bandwidth of the bandpass filter, which is the most important parameter in determining the number of points 
in parameter space, is strong enough that efficient detection is computationally intensive. Consequently, our principal 
conclusion is that, although the symmetry described in section || simplifies the problem from a formal point of view, 
it does so by strongly increasing the level of noise, principally by mixing components of noise of different frequencies. 
We have presented an algorithm that minimizes this effect, but have found that for a realistic broad-band search, it 
did not perform better than the stack-slide technique at the same computational power. 

Most detection algorithms for periodic sources, and in particular the stack-slide technique, probably don't have a 
strong scaling of their performances with the bandwidth of the search, essentialy because making use of the reduced 
bandwidth adds complexity by forcing one to consider the frequency as an explicit parameter of the search. It 
is therefore possible that our algorithm compares more advantageously to these other algorithms in narrow-band 
searches, which can be considered as more natural problems for it, because it necessarily involves a search over the 
source frequency. 
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